package optimizers.cmaes.dnacomputing;
import java.util.Random;

import model.Constants;
import model.SimpleSequence;

import visualization.GraphDisplayer;

import optimizers.cmaes.CMAEvolutionStrategy;
import optimizers.cmaes.fitness.IObjectiveFunction;
import optimizers.commons.GAIndiv;
import optimizers.tools.FileManager;
import optimizers.tools.Misc;
import fitness.BistableFitness;
import fitness.Fitness;
import gaforbistability.GAIndiv;

/*
 * 2013-02-07: handle lower bound value (replace x[index]<0?x=0:x=x[index] by x[index]*x[index] )
 */


/** DNA Computing interface for fitness evaluation 
 */
class DNACompFitness1 implements IObjectiveFunction { // meaning implements methods valueOf and isFeasible

	final static boolean COMPLEX_FITNESS = false;
	
	public static FileManager myLogFile;
	
	int numberOfSimpleSequences = -1;
	public int currentGen = 0; // not used for the moment
	public int currentInd = 0; // not used for the moment
	
	int evaluation = 0;
	
	double bestEverValue = Double.NEGATIVE_INFINITY;
	
	public DNACompFitness1(int __numberOfSimpleSequences) 
	{
		this.numberOfSimpleSequences = __numberOfSimpleSequences;
	}

	public void setInfo(int __currentGen, int __currentInd) // NOT USED FOR THE MOMENT
	{
		this.currentGen = __currentGen;
		this.currentInd = __currentInd;
	}
	
	// ***
	
	public double convertValue( double x )
	{
		double res;
		
		//res = x*x; // the advised method. HOWEVER this turns out to be computationally very extensive because of non-sparse matrix.
		//res = x<0?0:x ; // violent bound at zero, not classical, but in this particular case this may help as sparse matrix becomes easy to set. 
		
		if ( x > 0.5 ) // \_/, centered on 0.
			res = x - 0.5;
		else
			if ( x < -0.5 )
				res = -x - 0.5;
			else
				res = 0;
		
		return res;
	}
	
	// ***
	
	public double valueOf (double[] x) {

		double res = 0;
		
		double[][][] myTemplateValues = new double[numberOfSimpleSequences+1][numberOfSimpleSequences][numberOfSimpleSequences];
		double[][] mySequenceValues = new double[numberOfSimpleSequences+1][numberOfSimpleSequences];
		double[] mySeqK = new double[numberOfSimpleSequences];
		
		double exo;
		
		int index = 0;
		
		double normalizeValue;
		
		// [!!!]:
		//   - in the following, genome values are rescaled to the domain required + lower bound (x>0 mandatory for every dimension) is addressed.
		//   - "a lower bound a for a variable x is to use a + x2 instead of x in the objective function call" (cf. cmaes site)
		
		normalizeValue = 200;
		
		for (int i = 0 ; i < numberOfSimpleSequences+1 ; i++ ) // concentration des templates
			for (int j = 0 ; j < numberOfSimpleSequences ; j++ ) 
				for (int k = 0 ; k < numberOfSimpleSequences ; k++ )
				{
					myTemplateValues[i][j][k] = convertValue(x[index]) * normalizeValue;  // between 0 and 200
					index++;
				}
		
		normalizeValue = 100;
		
		for (int i = 0 ; i < numberOfSimpleSequences+1 ; i++ ) // concentration initiale : concentration sequence simple (numberofSimpleSeq) + sequence de tous les inhibiteurs (numberofSimpleSeq*numberofSimpleSeq)
			for (int j = 0 ; j < numberOfSimpleSequences ; j++ )
			{
				mySequenceValues[i][j] = convertValue(x[index]) * normalizeValue;  // between 0 and 100
				index++;
			}

		normalizeValue = 2e1/Constants.Kduplex - 3e-2/Constants.Kduplex;
		double offsetValue = 3e-2/Constants.Kduplex;
		
		for (int i = 0 ; i < numberOfSimpleSequences ; i++ ) // seqK : parametre thermodynamiaue
		{
			double value = convertValue(x[index]) * normalizeValue + offsetValue; // between 3e-2/Constants.Kduplex and 2e1/Constants.Kduplex
			
			if ( value < Constants.simpleKmin ) 
			{
				mySeqK[i] = Constants.simpleKmin;
			}
			else 
				if ( value > Constants.simpleKmax )
				{
					mySeqK[i] = Constants.simpleKmax;
				}
				else
					mySeqK[i] = value;
			
			index++;
		}

		normalizeValue = 1.0;
		exo = convertValue(x[index]) * normalizeValue; // between 0 and 1 // concentration enzyme exo-nuclease

		GAIndiv myIndividual = new GAIndiv(numberOfSimpleSequences, new Random(), myTemplateValues, mySequenceValues, mySeqK, exo, currentGen, currentInd); 
		
		myIndividual.OligoModelComplex.saturableExonuclease = false;
		myIndividual.OligoModelComplex.saturableNick = false;
		myIndividual.OligoModelComplex.saturablePoly = false;
		
		myIndividual.OligoModelComplex.seqK = mySequenceValues;
		myIndividual.OligoModelComplex.template = myTemplateValues;
		myIndividual.OligoModelComplex.exo = exo;

		//GraphDisplayer displayer = new GraphDisplayer("");
		
		//displayer.addGraph(myIndividual.printVisualModel(), myIndividual.displayTimeSereisWithInput(), null);
		//displayer.display();
		//System.out.println(myIndividual.evaluateFitness());
		
		//System.out.println(myIndividual.OligoModelComplex.printModeltoString(numberOfSimpleSequences, null));
		
		// [!!!] CMA-ES minimizes F ! (<==> maximize -F)

		Fitness myFitness = new BistableFitness();
		if ( COMPLEX_FITNESS == true )
			res = (myFitness.evaluateFitness(myIndividual)).finalResult(); // complex, computer-intensive fitness
		else
			res = myFitness.evaluateFitnessSimple(myIndividual).finalResult(); // Simple fitness
			
		// COMPUTE&APPLY PENALIZATION WRT SPARSITY
		double sparsity = 0;
		for ( int i = 0 ; i < x.length ; i++ )
			if ( x[i] > -0.5 && x[i] < 0.5 )
				sparsity = sparsity + 1;
		sparsity = sparsity / (double)x.length;
		res = res * sparsity;
		
		// monitor and display
		
		if ( res > this.bestEverValue )
			this.bestEverValue = res;
		
		System.out.println(this.evaluation + "," + this.currentGen + "," + this.currentInd  + "," + res + "," + this.bestEverValue);
		myLogFile.write(this.evaluation + "," + this.currentGen + "," + this.currentInd  + "," + res + "," + this.bestEverValue + "\n");

		this.evaluation++;
		
		return -res; // MIMIMIZE
	}
	
	// *** 
	
	public boolean isFeasible(double[] x) {return true; } // entire R^n is feasible
}

public class CMA_DNAComp_v1 {
	public static void main(String[] args) {
		
		int numberOfSimpleSequences = 4; //Matrix size is nxnx(n+1) -- value in [4...6]
		

		IObjectiveFunction fitfun = new DNACompFitness1(numberOfSimpleSequences);

		// new a CMA-ES and set some initial values
		CMAEvolutionStrategy cma = new CMAEvolutionStrategy();
		
		//cma.readProperties(); // read options, see file CMAEvolutionStrategy.properties
		cma.readProperties("src/optimizers/config/CMAES-DNA-20130207.properties");
		
		cma.setDimension( numberOfSimpleSequences*numberOfSimpleSequences*(numberOfSimpleSequences+1) + (numberOfSimpleSequences+1)*numberOfSimpleSequences + numberOfSimpleSequences + 1); // overwrite some loaded properties

		cma.setInitialX(0.0); // in each dimension, also setTypicalX can be used
		cma.setInitialStandardDeviation(1.5); // 0.15 // also a mandatory setting 
//		cma.options.stopFitness = 1e-14;       // optional setting

		// initialize cma and get fitness array to fill in later
		double[] fitness = cma.init();  // new double[cma.parameters.getPopulationSize()];

		String timeStamp = Misc.getCurrentTimeAsCompactString();

		// initial output to files
		//cma.writeToDefaultFilesHeaders(0); // 0 == overwrites old files
		cma.writeToDefaultFilesHeaders("logs-optimizers/log_"+timeStamp+"-",0);
		
		cma.helloWorld();
		
		// log info
		((DNACompFitness1)fitfun).currentGen = 0;
		((DNACompFitness1)fitfun).currentInd = 0;
		
		System.out.println("Time stamp: " + timeStamp);
		System.out.println("# evaluation,currentGen,currentInd,currentFitess,bestFitnessEver");
		
		((DNACompFitness1)fitfun).myLogFile = new FileManager("logs-optimizers/log-cmaes-"+timeStamp+".txt");
		
		((DNACompFitness1)fitfun).myLogFile.write("# evaluation,currentGen,currentInd,currentFitess,bestFitnessEver\n");


		// iteration loop
		while(cma.stopConditions.getNumber() == 0) {

            // --- core iteration step ---
			double[][] pop = cma.samplePopulation(); // get a new population of solutions
			for(int i = 0; i < pop.length; ++i) {    // for each candidate solution i
            	// a simple way to handle constraints that define a convex feasible domain  
            	// (like box constraints, i.e. variable boundaries) via "blind re-sampling" 
            	                                       // assumes that the feasible domain is convex, the optimum is  
				while (!fitfun.isFeasible(pop[i]))     //   not located on (or very close to) the domain boundary,  
					pop[i] = cma.resampleSingle(i);    //   initialX is feasible and initialStandardDeviations are  
                                                       //   sufficiently small to prevent quasi-infinite looping here
				//System.out.print(".");

				// compute fitness/objective value	
				fitness[i] = fitfun.valueOf(pop[i]); // fitfun.valueOf() is to be minimized
				
				((DNACompFitness1)fitfun).currentInd++;
			}
			
			//System.out.print("\n");
			
			cma.updateDistribution(fitness);         // pass fitness array to update search distribution
            // --- end core iteration step ---

			// output to files and console 
			cma.writeToDefaultFiles();
			int outmod = 150;
			if (cma.getCountIter() % (15*outmod) == 1)
				cma.printlnAnnotation(); // might write file as well
			if (cma.getCountIter() % outmod == 1)
				cma.println(); 
			
			((DNACompFitness1)fitfun).currentGen++;
			((DNACompFitness1)fitfun).currentInd = 0;
		}
		// evaluate mean value as it is the best estimator for the optimum
		cma.setFitnessOfMeanX(fitfun.valueOf(cma.getMeanX())); // updates the best ever solution 

		// final output
		cma.writeToDefaultFiles(1);
		cma.println();
		cma.println("Terminated due to");
		for (String s : cma.stopConditions.getMessages())
			cma.println("  " + s);
		cma.println("best function value " + cma.getBestFunctionValue() 
				+ " at evaluation " + cma.getBestEvaluationNumber());
			
		// we might return cma.getBestSolution() or cma.getBestX()

		((DNACompFitness1)fitfun).myLogFile.close();
		
	} // main  
} // class



